WhereWulff: A Semiautonomous Workflow for Systematic Catalyst Surface Reactivity under Reaction Conditions

This paper introduces WhereWulff, a semiautonomous workflow for modeling the reactivity of catalyst surfaces. The workflow begins with a bulk optimization task that takes an initial bulk structure and returns the optimized bulk geometry and magnetic state, including stability under reaction conditions. The stable bulk structure is the input to a surface chemistry task that enumerates surfaces up to a user-specified maximum Miller index, computes relaxed surface energies for those surfaces, and then prioritizes those for subsequent adsorption energy calculations based on their contribution to the Wulff construction shape. The workflow handles computational resource constraints such as limited wall-time as well as automated job submission and analysis. We illustrate the workflow for oxygen evolution reaction (OER) intermediates on two double perovskites. WhereWulff nearly halved the number of Density Functional Theory (DFT) calculations from ∼240 to ∼132 by prioritizing terminations, up to a maximum Miller index of 1, based on surface stability. Additionally, it automatically handled the 180 additional resubmission jobs required to successfully converge 120+ atoms systems under a 48-h wall-time cluster constraint. There are four main use cases that we envision for WhereWulff: (1) as a first-principles source of truth to validate and update a closed-loop self-sustaining materials discovery pipeline, (2) as a data generation tool, (3) as an educational tool, allowing users (e.g., experimentalists) unfamiliar with OER modeling to probe materials they might be interested in before doing further in-domain analyses, (4) and finally, as a starting point for users to extend with reactions other than the OER, as part of a collaborative software community.

One can compute the surface energy of a symmetric slab: distributing it between the two surfaces, each having area A, where the sum is over the chemical potentials µ i and the number of atoms N i of each of the species that make up the composition of the slab, as shown in Eq. (1).
Assuming that the bulk is in equilibrium with the slab and by defining a reference, we can conveniently re-write Eq.
(1) to be as a function of the bulk energy per formula unit, as shown in Eq.
(2), where the last sum represents the free energy excess, x i the number of atoms per bulk formula and N ref is the reference specie that is picked.
More concretely, the derivation of Eq.
(2) for a slab consisting of N Ba , N T i , N O and a bulk with formula BaTiO 4 : • From first principles we can express the Gibbs energy of the system as the sum of the interfacial energy and the cost of introducing/removing elemental species i at some temperature (T) and pressure (p) • Rearranging Eq. (3) for the surface energy, γ, and assuming symmetric surfaces, we get Eq. (4) S-3 • If we assume that the bulk is in equilibrium with the slab, we can write Eq. (5) • We define the reference to be in relation to Ba per Eq. (6) • We then introduce the bulk energy per formula unit into Eq. (4) by substituting Eq. (6), to get Eq. (7) • Eq. (7) reduces to the well-known surface energy Eq. (8) when the stoichiometry between the bulk and the interface is maintained The slab model is made up of a 2D surface and a corresponding oriented bulk, since this has been shown to most efficiently converge the surface energy calculations. S1

S2 Surface Coverage and OER
The coverage effects of reaction intermediates (OH*, O*) may significantly impact the local environment of the active site, resulting in changes in the adsorption strength of the OER reaction intermediates. Consequently, it is important to determine the most stable surface coverage at given conditions of applied potential and pH to more reasonably describe the OER activity. To determine the most stable surface coverage, Surface Pourbaix diagrams, S2 were constructed at three extreme coverages, clean, OH*, and O* terminated.
OER catalytic activities of the different surface structures were determined by the theoretical S-4 overpotential (η OER ) and the potential-determining step (PDS) by assuming an associative reaction mechanism with O*, OH* and OOH* as reaction intermediates. The four proton-coupled electron transfers (PCET) reactions under acidic conditions are: The computational hydrogen electrode (CHE) was used to express the chemical potential of The theoretical thermodynamic OER overpotential, which is a measure of the activity of a catalyst, is then defined from Gibbs free energies of reactions (9)-(12): The step with the largest value in Eq. (13) is referred to as the potential-determining step (PDS).
It is important to note that the overpotential should not be compared directly with a measured S-5 overpotential, since the measured overpotential depends on the current density. Figure S1: Standard convergence tests on bulk BaSnTi 2 O 6 illustrating the trade-off between accuracy and computational time. Based on these results, we infer that using an ENCUT setting of 500 eV and a kpoint grid of 5×5×5 would be a good balance between accuracy and runtime. Our custom DFT settings, implemented as a child class of Pymatgen's MVLSlabSet, obey these thresholds.

S3 Convergence Tests and Settings
S-6  (10). These results corroborate that it is infeasible for us to carry out the slab optimizations at those levels of thickness and under such large numbers of degrees of freedom, which are shown to trigger multiple re-submissions ( Figure S2).